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^ ; Abstract 

We describe discretisations of the shallow water equations on the sphere using the frame- 
work of finite element exterior calculus, which are extensions of the mimetic finite difference 
framework presented in (Ringler, Thuburn, Klemp, and Skamarock, Journal of Computa- 
tional Physics (2010)). We present two formulations, a "primal" formulation in which the 
<^ ■ finite element spaces are defined on a single mesh, and a "primal-dual" formulation in which 
a finite element spaces on a dual mesh are also used. Both formulations have velocity and 
layer depth as prognostic variables, but the exterior calculus framework leads to a conserved 
diagnostic potential vorticity. In both formulations we show how to construct discretisations 
that have mass-consistent (constant potential vorticity stays constant), stable and oscillation 
free potential vorticity advection. 
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1. Introduction 



In a recent paper on horizontal grids for global weather and climate models, [ST12I ] listed 
a number of desirable properties that a numerical discretisation should have, which can be 
paraphrased as accurate representation of geostrophic adjustment, mass conservation, curl- 
free pressure gradient, energy conserving pressure terms, energy conserving Coriolis term, 
steady geostrophic modes, and absence/control of spurious modes. Of this list as presented 
here, the first property could be said to relate to the stability and accuracy of the discrete 
Laplacian formed from divergence and gradient operators, whilst the next five all relate to 
mimetic properties [i.e. the numerical discretisations exactly represent differential calculus 
identities such as V x V - 0), and the last prop erty relates to the kernels of the various discre- 
tised operators (see |LSRH05l . ILRP07I . |LP08| and related papers by Le Roux and coworkers 



for extended discussion of these issues in the context of finite element methods). In the 
context of the rotating shallow-water equations on the sphere, which represent the standard 
nonlinear framework for investigating horizontal grids for global models, the C-grid stagger- 
ing on the latitude-longitude grid combined with an appropriate choice of reconstruction of 
the Coriolis term provides all of these properties, but leaves us with a grid system with a 
polar singularity. This, together with a need for models with variable resolution, has started 
a quest for alternative grids and discretisations that satisfy these properties. 
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The extension of the C-grid to triangular meshes (and the finite element analogue, the 
RTO-PO discretisation) satisfies the firs t six pr operties and has been popular in both atmo- 



sphere and ocean applications f |WC98l . |BR05| ). however it is now well understood that the 



triangular C-grid supports spurious inertia-gravity mode branches because of the decreased 
ratio of vel ocity degrees of freedom (DOFs) to pressure DOFs relative to quadrilaterals (from 
2:1 to 3:2) |DanlO . Gasll |. More recently, a Coriolis reconstruction for the hexagonal C-grid 
was derived in Thu08l ] that provides the mimetic properties described above, and this was 
extended to arbitrary orthogonal C-grids (grids in whi ch dual g rid edges that join pressure 
points intersect the primal grid edges orthogonally) in |TRSK09 ]. The hexagonal C-grid has 
an increased ratio of velocity DOFs to pressure DOFs (from 2:1 to 3:1), and so does not 
support spurious inertia-gravity mode branches, but does have a branch of spurious Rossby 
modes. This reconstruction can be used to construct energy and enstrophy conserving C-grid 
discre tisations for the nonlinear rotating shallow-water equations using the vector invariant 



form [RTKSlOj . in which mimetic properties are used to produce a velocity-pressure formu- 
lation in which the diagnosed potential vorticity is locally conserved, is consistent with the 
mass conservation (i.e. constant potential vorticity stays constant in the unforced case), and 
can be controlled such that it remains bounded. 

Two directions remain outstanding from this approach, namely the relaxation of the or- 
thogonality requirement which constrains cubed sphere grids so tha t grid resolution increases 
much more quickly in the corners that at the middle of the faces PL07I ] . and the construc- 
tion of higher-order operators to av oid gri d imprinting. Two recent papers by the authors 
attempted to address t his issues. In |TC12j | a framework was set up to generalise the mimetic 
approach of RTKSld| to non-orthogonal grids, b ut the method of constructing sufficiently 
high-order operators was not clear. Meanwhile, in CS12I ]. it was show n that mi xed finite ele- 
ment methods in the framework of finite element exterior calculus (see AFW06| for a review) 
provide the first six properties listed above, plus sufficient flexibility to adjust the ratio of 
velocity DOFs to pressure DOFs to 2:1 to avoid spurious mode branches. The BDFM1 space 
on triangles and the RTk heirarchy of spaces on quadrilaterals were advocated as examples 
of spaces that satisfy that ratio. However, in that paper it was not clear how the extension 
to nonlinear shallow-water equations would be made. In this paper we address both of these 
open questions by describing a finite element exterior calculus framework for the shallow- 
water equations, which enables us to write the equations in a very compact form and reveals 
the underlying structure behind the mimetic properties. We shall discuss two formulations, 
a primal grid formulation in which potential vorticity is represented on a continuous finite 
element space, and a pr imal-dual grid fo rmulation that makes use of the discrete Hodge star 
operator introduced in HipOlal . iHipOlbl ] , in which potential vorticity is represented on a dis- 
continuous finite element space. In the latter case, discontinuous Galerkin or finite volume 
methods can be used for locally conservative, bounded, mass-consistent potential vorticity ad- 
vection, whilst in the former case we show that streamline-upwind Petrov- Galerkin methods 
with discontinuity-capturing schemes can be be incorporated into the framework to provide 
conservative, high-order, stable, non-oscillatory advection of potential vorticity. 

The rest of this paper is structured as follows. In Section [2] we provide a "hands-on" 
introduction to the calculus of differential forms, then write the rotating shallow water equa- 
tions in differential form notation. In Section 13.21 we describe our primal grid finite element 
exterior calculus formulation of the shallow water equations, and in Section @] we describe 



2 



our primal/dual grid formulation. Finally, in Section we provide a summary and outlook. 



2. Differential forms on manifolds 

In this section, we introduce the required elements from the language of differential 
forms, i n an informal manner. For more rigorous definitions, the reader is referred to 



MR991 IAFW06I . lHolll| . We then combine these elements to write the rotating shallow- 



water equations on the sphere in differential form notation. 
2.1. Differential form preliminaries 

Solution domain. We shall consider the case in which the solution domain Q is a closed 
compact two-dimensional surface embedded in three dimensional space. In applications the 
main surfaces of interest are the surface of the sphere, or a rectangle in the x-y plane with 
periodic boundary conditions in both Cartesian directions. For brevity of notation we do 
not consider domains with boundaries; this avoids the need to include boundary terms when 
integrating by parts, although they can easily be included. We shall expand all quantities in 
three-dimensional Cartesian coordinates, and shall use the machinery of contraction operators 
(defined later) to restrict the equations to the two-dimensional surface Q. Hence, we begin 
by developing differential calculus in R 3 . 

Vector fields. Informally, a vector field u on R 3 is a mapping u : R 3 -*■ R 3 , i.e. a vector field 
assigns arbitrary 3-vectors to each point in R 3 . We denote 3C(R 3 ) as the space of vector fields 
tangent to Lagrangian flow lines on R 3 . Physically, velocity fields are represented as vector 
fields. In due course we shall also require vector fields on Q, which are only defined for points 
x € Q, and the assigned vectors are constrained to be tangential to Q. We denote as 
the space of vector fields on Q. 

Differential forms. In this paper we shall make use of four types of differential forms: 0- 
forms (which are simply scalar- valued functions), 1-forms, 2-forms, and 3-forms. In general, 
1-forms are used to compute line integrals, 2-forms are used to compute surface integrals, and 
3-forms are used to compute volume integrals. We shall denote A k as the space of A;-forms. 

Cotangent vectors on R 3 . To define differential forms, we first need to define cotangent vectors 
which are the dual objects to vectors, i.e. mappings from tangent vectors to R. In R 3 , 
cotangent vectors have the basis {dx l } 3 =1 which is dual to the Cartesian basis {x 1 }™^ for 
tangent vectors, so the general form for cotangent vectors is 



m 



^rriidx 1 := m ■ dx, m= (m 1 ,m2,m 3 ) T , dx = (dx l ,dx 2 ,dx 3 ) T . 



Cotangent vectors m are defined by their result when the inner product is taken with tangent 
vectors v, i. e. (m ■ dx,v) = m T v. 
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1-forms onR 3 . A differential 1-form uj assigns a cotangent vector to each point x e M 3 . The 
inner product (u, u) of a 1-form uj and a vector field u is therefore a scalar function (0-form) 
defined on M 3 . On M 3 , it is convenient to expand 1-forms in the standard basis for three- 
dimensional cotangent vectors, i.e., (dx 1 , dx 2 , dx 3 ). Informally, a 1-form can be integrated 
along a one-dimensional curve CcR 3 using the usual definition of line integration 

Wedge product and volume forms. The wedge product can be used to construct 2-forms and 
3-forms from 1-forms. The wedge product of a scalar function (0-form) / with a A;-form u is 
simply the arithmetic product: 

/ auj = fu. 

In general, the wedge product of an n-form and an m-form produces an (n+m)-form. The 
wedge product satisfies the following conditions: 

1. Bilinearity: 

(aa.\ + ba.2) auj = a(«i auj) + b(a 2 au), a a (aa>i + buj 2 ) = a(a a uj\) + b(a a uj 2 ), 

2. Anticommutativity: 

UA7 = (-1) H 7AW, 

for a /c-form ui and an /-form 7, and 

3. Associativity: 

(u) A 7) A K = OJ A (7 A k). 

From these properties it may be deduced that the wedge product of two arbitrary 1-forms cu, 
7 takes the form 

uj a 7 = ct\ dx 1 a dx 2 + a 2 dx 1 a dx 3 + a 3 dx 2 a dx 3 , 

for some scalar functions a±, a 2 , «3, and hence this is the general form for 2-forms in Cartesian 
coordinates. Similarly it may be deduced that the wedge product of three arbitrary 1-forms 
uj, 7 and k takes the form 

uj a 7 a k, = f(x) dx 1 a dx 2 a dx 3 , 

for some scalar function f(x), and hence this is the general structure of 3-forms. We define 
the "volume form" 

dV = dx 1 a dx 2 a dx 3 . 

A general 3-form f(x)dV may be integrated over three-dimensional submanifolds M e M. 3 
using the usual definition of triple integration: 

f fdV= f fdx 1 dx 2 dx 3 . 

JM JM 
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Contraction with vector fields. Having denned the volume form, we can use contraction with 
vector fields to construct 1-forms and 2-forms on the surface Q. In general, contraction of a 
vector field u with a /c-form u results in a (k-l)-form, denoted u j oj. The contraction of a 
vector field u with a 1-form w is simply the inner product 

U J UJ - (u, cu) - Y] U l U)i. 

i 

The contraction is linear, i.e. u j (a(x)ui + b(x)ui2) - a{x)u J uj\ + b(x)u j for two scalar 
functions a(x) and b(x), and two A;- forms Wi and Hence, we may define the contraction 
with 2-forms for the basis elements dx 1 a dx 2 , dx 1 a dx 3 , and dx 2 a dx 3 , with 

u j d x' 1 a d x 3 = Ui d x 3 - uj d x l , 1 < i + j < 3. 

Similarly, the contraction with 3-forms can be defined from 

u j d x 1 a d x 2 a d x 3 = U\ d x 2 a d x 3 - u 2 d x 1 a d x 3 + % d x 1 a d x 2 . 

Surface form on Q. The surface form dS is used for surface integrals over two-dimensional 
submanifolds M of f2; it is defined on Q &s d S = k -i dV where k is the unit oriented vector 
field normal to Q. The general form of 2-forms on Q is then 

uj = p(x) d S, 

for scalar functions p(x) defined on Q. A 2-form on Q can be integrated over a two dimen- 
sional submanifold M of Q using the usual definition of surface integration: 

f u= [ p(x)dS. 

Identification of vector fields on Q with 1-forms. In this framework we shall make use of two 
different identifications of vector fields on Q with 1-forms. The first is written 

u h+ u- dx n = Y,Uidx l n , 

i 

where 

dxn = dx - k(k ■ dx), 

i.e. the projection of dee into the dual space to the tangent plane on Q at x. This identifi- 
cation is used to compute circulation integrals 

/ u- dxQ= / u • dec, 
Jc Jc 

along curves C c fi, and hence is associated with the curl operator on Q. The second 
identification is written using the contraction with dS, 

u h+ u j d S, 

and is used to compute flux integrals 

/ u-idS= / u-kxdx 
Jc Jc 

across curves Ccfi, and hence is associated with the divergence operator on Q. 

We shall denote A°(f2) as the space of 0-forms (scalar functions) on Q, A 1 (fi) as the space 
of 1-forms on Q, which can be written in both of the identifications with vectors fields given 
above, and A 2 (Q) as the space of 2-forms on Q written in the form u = p(x) dS where p is 
a scalar function. 
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Differential operator. The differential operator d neatly encodes all of the vector calculus 
differential operators e.g., div, grad, curl etc. In general, the differential operator d maps 
/c-forms to k + 1-forms, and satisfies: 

1. For scalar functions /, d / = £i §£r dx l = V/ • dx. 

2. Product rule: d(u a 7) = (duo) a 7 + uj a (d7). 

3. Closure: <=> d(daj) = d 2 a; = 0. 

We shall see some examples very shortly, however, first we shall write the general form of 
Stoke's theorem for d: 



f du= f 

JM Jd 



CO, 

dM 



where M is a /c-dimensional submanifold of M 3 , uj is a (A;-l)-form (and hence da; is a fc-form), 
and dM is the (k - l)-dimensional submanifold of M. 3 corresponding to the boundary of M. 
Combining Stoke's theorem with the product rule provides the integration by parts formula 



/ (da;)A7= / uj a 7 - / 0JA(d7). 

JM JdM JM 



Differential operators on Q. Standard vector calculus differential operators on scalar func- 
tions / and vectors fields u defined on Q are obtained from the two identifications of vector 
fields with 1-forms: 

df = Wf-dxn, 

df = (fcxV/)jdS:=v7jd5, 
d(u,-da^) = k ■ V x udS := V 1 ■ udS, 
d(uJdS) = V-udS, 

where V, V 1 , V 1 - and V- are differential operators defined intrinsically on the two-dimensional 
surface Q. The closure property d 2 = then leads to the following vector identities for the 
two identifications of vector fields with 1-forms: 

= d 2 / = d(V/d : i ;n ) = V 1 -V/d 1 S, (1) 
= d 2 / = d(vV JdS) = V- V 1 /^. (2) 

These identities are crucial for geophysical applications since they dictate the scale separation 
between slow and fast dynamics. 

Hodge star on Q. The Hodge star operator * on Q maps from /c-forms to (2-&)-forms. It is 
used in this paper to define the L 2 -inner product between two /c-forms u and 7 by 



<w,7>= f 
Jn 



U A *7, 

u 

and is also used to define the Coriolis term. The Hodge star is linear (i.e., *(a(x)u + b(x)-y) = 
a(x) * u + b(x) * 7 for scalar functions a, b and /c-forms uj, 7), and so we may define it for our 
chosen Cartesian basis, where the Hodge star satisfies 

*/ = fdS, 

*fdS = f, 

*(u-dxfi) = u j dS = k x u ■ dxn, 

*u-idS = -(u ■ dxn) = k x u j dS. 
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From the presence of kx in these formulas it becomes clear why the Hodge star might be 
relevant to the Coriolis term. Note that ** = Id for 0- and 2-forms, and ** = - Id for 1-forms. 

2.2. Rotating shallow-water equations in differential form notation 

We have now established enough notation to write the rotating shallow-water equations 
on Q in differential form notation. We start from the following form of the rotating shallow- 
water equations: 

^-u+(( + f)u' + V(g(D-b) + ±\u\ 2 ) = 0, 

1d + V-(uD) = 0, 

where u is the velocity, C = fc-Vxn = V i -uis the vorticity, D is the layer depth, b is the 
height of the bottom surface, g is the acceleration due to gravity, / is the Coriolis parameter 
and 1 = kx. 

Using the notation that we have described above, we can immediately rewrite these 
equations as 

-u-dx n + *((( + f)u-dx n ) + d(g(D-b) + -\u\ 2 ] = 0, (3) 

^DdS + d(u jDdS) = 0, (4) 
d(«-dx n ) = CdS. (5) 

Applying d to Equation ([3]) and making use of d 2 = and the definition of Hodge star 
immediately gives the vorticity equation 

^CdS + d(uj(C + /)dS) = 0, 

which is in the same flux form as the mass equation (Equation (jll)). The potential vorticity 
(PV) q is defined from 

(C + f)dS = qDdS, 
and hence we obtain the law of conservation of potential vorticity 

^(qDdS) +d(u ^qDdS) = 0. (6) 
Note that if q is constant then 

/ 



^DdS + d(uJqDdS) 



dq 



\ .0 / 

which means that q remains constant. This is what is meant by consistency with Equation 

Our goal is to design a framework for finite element discretisations that has u and D as 
the prognostic variables, yet preserves the conservation law structure of equations (@J and 
([6]). Furthermore, we shall show how stabilisations for these conservation laws (which are 
required for meteorological applications) can be incorporated into this framework. 
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3. Finite element exterior calculus formulation 



3.1. Finite element spaces 

The fundamental idea of finite element exterior calculus applied to the rotating shallow- 
water equations is to choose finite element spaces for u, £ and D such that the operator d 
maps from one expansion to another, so that the vector calculus identities ([1]) and ([2]) still 
hold. The difficulty is that continuous derivatives in the normal direction across element 
boundaries are required to compute d(u -idS), whilst continuous derivatives in the tangent 
direction across element boundaries are required to compute d(u ■ dx^i). On a single grid, 
we cannot have both, and thus must choose to construct finite element expansions such that 
only one of ([!]) or (T5]) hold in the strong form, and the other will hold in the weak form 
having integrating by parts. This amounts to choosing one of D and ( to have a continuous 
finite element expansion and the other to have a discontinuous expansion. A discontinuous 
expansion allows for discontinuous Galerkin methods which are locally conservative and allow 
for monotonic advection schemes, and in meteorogical applications it is more important that 
these schemes are available for D than £, so we choose to hold ([[]) in the strong form. Later, 
in setting up the primal-dual grid formulation, we shall introduce a dual grid for which (T21) 
holds in the strong form, consistently with the weak form on the primal grid. 

Finite element differential form spaces. Having made the choice to hold (CQ) in strong form, 
we need to choose finite element spaces for £, u, and D, denoted E, S and V respectively. 
This choice defines equivalent subspaces A k (Q) c A k (Q), k = 1,2,3, given by 

A°(fi) = E, 

k\n) = {uJdS:ueS}, 
A 2 (Q) = {DdS-.DeV}, 

We require that d maps from A 1 (f2) into A 2 (f2), and that d maps from A°(f2) into A 1 (r2), 
which implies that E is a continuous finite element space, S is div- conforming (i.e. u e S 
has continuous normal components across element boundaries), and V is a discontinuous 
finite element space. Numerous examples of (E, S, V) satisfying these properties, including 
E = P(k), S = BDM(k), V = P{k) DG . 

Dual differential operator. We define 6 as the dual differential operator from A fc to A^ 1 that 
is dual to d, i.e. 

/ 7 a *5uj = - / d^j a *oj, V7 e A fc . 

We note that 5 2 u = 0, since 



/ 7 a *i 2 u = - d~f/\*5u 
Jn Jn 

/ d 2 7A*w 
Jn 

= V7eA°,u;eA 2 . 
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Discrete Helmholtz decomposition. As discussed in AFW06], if d maps from A (Q) onto 
the kernel of d in A 1 (f2), then there is discrete Helmholtz decomposition and any 1-form 
u -idS e A 1 (f2) can be written as 

u j dS = dip + 5(f) + P), 

where ip e A°(f2), e A 2 (f2), and f) e "H, where % c A 1 (fi) is the space of discrete harmonic 
1-forms given by 

^ = {f)eA 1 (^):d() = 0, 5t) = 0}, 

which has the same dimension as the space of continuous 1-forms on Q (which has dimension 
for the surface of a sphere). 

Construction of global finite element spaces by pullback. We construct the spaces A fc (f2), k = 
0,1,2, element by element, and then enforce continuity constraints, namely: 7 e A°(f2) is fully 
continuous, u J dS e A x (f2) has continuous normal components across element boundaries, 
and (ftdS e A 2 (f2) is allowed to be fully discontinuous. The restriction A fc (e) to element e of 
A k (Q) is defined through a pullback mapping e : A fc (e) -> A fc (e) from a reference element e 
(where integrals are computed) to each element e, the pullback 0* : A fc (e) -> A fc (e), is defined 
by 

/ / , 

~/M J(f>e(M) 

for all o> e A fc (e) and all integrable submanifolds Mee. Specifically: 

7 = 0* 7 = 7 o0, 7 eA°(e), 7 e A°(e) 

u J dS = (f)* e u j dS = it j fc j d0i a d0 2 a d0 3 , it J dS e A°(e), u J dS e A°(e), 

= ujd0, u-idSek 1 (e), 

DdS = (p* e DdS = DdS, 

where dS = dx a dy is the surface form on e with coordinates x, and (pi, i = 1,2,3 are the 
components of 0. In coordinates, the pullback of u j d S defines the Piola transformation 

(P*(U jdS) = 1 - 7 r--^UJdS, 

det(g)ft* 

and the pullback of D d S 1 defines the scaling transformation 

(/)*(DdS) = ^jrrrdS. 

The pullback operator 0* commutes with d, ie. d(p*u = (p*du, so it is sufficient to check 
that d maps from A fc (e) to A fe+1 (e), k = 1, 2 to guarantee that it maps from A fc (f2) to A k+1 (Q). 

3.2. Primal grid formulation 

Here we shall concentrate on the semi-discrete continuous time equations. 
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Discrete velocity equation. Since we have chosen to use the strong definition of d(u J dS), 
we need to transform equation ([3]) into a form where this choice of identification with vector 
fields is used. This is done by applying the Hodge star: 

^-u jdS + *(Q JdS) + *d*(g(D-b) + K)dS = 0, Q = (( + f)u, 

where KdS = \u\ 2 /2dS and we have replaced D - b + K with *(D - b + K) dS. Note that 
we shall need to be careful with our definition of the PV flux Q for consistency, this shall be 
discussed later. 

Since A 2 (f2) is a discontinuous space, we cannot apply d to DdS and must adopt the 
weak form. This is done by taking the L 2 inner product with a test 1-form w J d S using 
the Hodge star, and integrating by parts (with vanishing boundary term since there are no 
boundaries): 

4- f wjdSA*ujdS- f wjdSAQjdS- f d(wjdS)A* (g(D - b) + K) dS = 0, VwJdS e A 1 ^). 
dt Jn Jn Jn 

(7) 

To obtain the Galerkin finite element approximation of this equation we simply restrict 
u JdS 1 , w j dS and Q J dS to the finite element space A 1 (f2), and DdS and bdS to the 
finite element space A 2 (Q). 

Discrete mass equation. In deriving the equation for D, the integral must be performed over 
a single element e due to the discontinuity, following the discontinuous Galerkin approach. 
Taking the L 2 inner product of a test 2- form (j)d S e A 2 (f2) with equation (j4j) over one element 
e gives 

— f(f)dS A*(Dds) = - f(f)dS A*d(ujDdS). 

dt Je Je 

To obtain coupling between elements we integrate by parts to obtain 
4 [<pdSA*(Dds) = fd*((j)dS)Au^DdS- f *(<f)dS) au jDdS, \/<f)dS. (8) 

d t Je J e J de 

where de is the boundary of element e and D is chosen as the value of D on the upwind side, 
following the standard discontinous Galerkin approach. Again we obtain the finite element 
approximation to equation fl4]) upon restriction of <j)dS and DdS to A 2 (f2), and u J dS* 
to A 1 (fi). In the case of A 2 (f2) being piecewise constant functions, we may select 4> as the 
indicator function for element e, i.e., 



<Kx) 

and Equation ([8]) becomes 



1 if x 6 e, 
otherwise, 



4 [Dds = - [ ujDdS, 

dt Je Jde 



which is the usual finite volume scheme for D, in which case a higher-order upwind flux uD 
is required. 
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Discrete vorticity equation. Next we define the vorticity £. This must be done weakly by 
introducing a test function to Equation (E]) and integrating by parts, since d(u-dxn) cannot 
be evaluated in strong form. We obtain C <e A°(f2) from the solution of 

f 7A*C = " f d 7 A*nJd5, V 7 eA°(fi). 
Jn Jn 

Since d 7 e A x (f2) for arbitrary 7 € A 2 (f2), we may select w J dS* = -d 7 in Equation ([7]) to 
obtain 

— T 7 a*C = -T7 [ d 7 A*ujdS 
dt Jn dt Jn 

= Td 7 AQjd5+ f d 2 ^A(g(D-b) + K)dS, 
Jn Jsi^ 

=0 

= [d-fAQjdS, V 7 eA°(fi). 

This is the standard continuous finite element discretisation of the vorticity transport equa- 
tion. Conservation of vorticity is a direct consequence of this, upon choosing 7 = 1: 

— F(dS = - f d(l) aQ jdS = 0. 
dt Jn jf!^^ 

=0 

Discrete potential vorticity. Having diagnosed a vorticity, we can diagnose a potential vor- 
ticity q e A°(f2) from 

f 1 AqDdS= f 7 A(C + /)d5. 
Jn Jn 

Then, the discrete potential vorticity satisfies 

— fjAqDdS+ f d 7 aQ jdS = 0. (9) 
dt Jn Jn 

It remains to define Q so that the diagnosed potential vorticity advection is mass-co nsistent, 
stable and oscillation-free. We shall follow identical steps to those taken in RTKSlOj , namely 
we shall first diagnose a mass flux F from Equation (jHJ), and then shall insert Q = Fq in 
Equation (jHJ), before making appropriate modifications to stabilise the advection. 

Discrete mass flux. To compute the mass flux, the crucial observation is that since d maps 
from A 1 (f2) onto A 2 (f2), there must exist a discrete mass flux F J dS e A 1 (SI) such that 

^Dds = -d(F jdS). (10) 

It turns out that F can be calculated locally, i.e. independently for each element without 
solving a global elliptic problem. Substitution of equation ( TTOl) into equation (JE]) and selection 
of (j) d S = d S gives 

fd(FjdS)= [ FjdS= [ uDjdS, 

Je J de J de 



11 



which is satisfied if we require that 



fde 



>dSA*F jdS 



f (f)dS a*uD jdS = 0, 

Jde 



€A 2 



The remaining degrees of freedom for F can be fixed by 
[d*((f)dS)AF jdS = - fd*((f)dS)Au jDdS+ f *((f)dS) a u j DdS, V0eA 2 , 

Je Je Jde 

plus local g auge co nditions that do not alter du jdS. This is essentially the Fortin projection 
into A 1 (ft) |For77 . 

This calculation can be extended to the time-discretised case in which a slope limiter is 
applied before each timestep or Runge-Kutta stage CSOJJ- In each element E, slope limiters 
aim to preserve monotonicity by adjusting DdS in each element E in such a way that De 
is preserved, where 

J E DdS 



D t 



f E dS 



We write the action of the slope limiter on DdS as 

S(DdS) = DdS + d(F s jdS). 
Since the slope limiter preserves the mean value of DdS, we obtain 

f F s jdS = 0, 

JdE 

which means that we can set 

£<f)F s jdS, V0eA 2 (ft), 
for all edges e € dE, and solve a problem for F s in the interior. 

where F is the discrete mass 



Discrete vorticity flux. We now choose Q = Fq in Equation 
flux defined above, and we obtain 



_d 
dt 



fjqDdS- [ d 7 AFq^dS = 0, V 7 e A . 
Jn Jn 



Note that since 7 and q are continuous, and F has continuous normal components, we may 
integrate by parts in the discrete equations to obtain 

4- [ iqDdS+ [~fAd(Fq JdS) = 0, V7 e A , 
dt Jn Jn 

We can now verify the consistency property, namely that constant q remains constant, by 
choosing q = 1 in the above equation, obtaining 



/ 



Jn dt Jn 



7 A 



\ 



0_ 

dt 



DdS + d{F JdS) 



V 



=0 



0. 



V 7 e A c 



dq 
di 



0. 



After computing Q, we can insert it into the velocity equation and we obtain prognostic 
velocity evolution that is consistent with this diagnostic PV flux. 
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Stabilised vorticity flux. As previously discussed, Equation (ITT]) is the continuous finite ele- 
ment approximation to the potential vorticity transport equation and thus requires stabili- 
sation. This can be done by introducing diffusive fluxes into Q J d S of an appropriate form, 
provided that the consistency property is preserved. The most basic type of stabilisation 
is from Laplacian diffusion that results from a down-gradient PV flux. Equation (ITT]) then 
becomes 

4- f jqDdS- f d 7 A (F jdS + *ndq) = 0, V 7 € A , 
dt Jn Jn 

where k is the diffusivity, and we obtain Q -i dS = qF J dS - *ndq. Note that this is 
mass-consistent since the diffusive flux vanishes for constant q. 

A m ore sop histicated approach is provided by the streamline-upwind Petrov-Galerkin 
method 



BH82 



which replaces the test function 7 with 7 + ru ■ d 7 , where r is a suitable 
scalar coefficient that depends on u and the mesh size h; this effectively biases the test 
function in the upwind direction. Substitution into Equation (ITT]) gives 



d_ 

dt 



fjqD 
Jn 



dS 



f d"f a F . 
Jn 



dS+ fd-fArujl^-(qDdS) + qFjds)=0, V 7 e A c 
Jn \ot 1 



and we obtain 



Q JdS = qF Jd5 + ™ J U^{qDdS) + qF Jds) 



As shown in [BH82[, the Petrov-Galerkin method is stable and converges at the optimum rate 
since it is a residual-based modification of the equations. It is still possible for overshoots and 
undershoots to appear near t o unres olved fronts, but these can be controlled by including a 
discontinuity capturing term HM86| . which consists of an additional Laplacian diffusion with 



spatially varying k which is non-zero only in regions where the solution is poorly resolved, 
in which case the flux becomes 



Q 



dS = qF jdS + TU^(-^(qDdS) + qFjds)-*Kdq. 



4. Primal-dual grid finite element formulation 

In this section we provide an alternative formulation that makes use of a second set of 
spaces defined on a dual grid based on the u ■ d identification of vector fields with 1-forms. 
The idea is that when we want to apply V 1 and V- operators strongly we use the primal grid 
spaces as defined in the previous section, and when we want to apply V and V 1 - strongly 
we use the dual grid spaces. This requires defining mappings between the primal and dual 
spaces which are defined via the Hodge star operator. 

We start by selecting a set of finite element differential form spaces A^ c A fc and A^(Q) c 
A fc on the primal grid and the dual grid respectively, satisfying 

Ag(fl) Aj(n) ki(n) 
AS(n) -?-> Aj(n) A3(n). 

We shall calculate with mass -DdS" in A^(fi), vorticity (dS in A^(f2), which are both dis- 
continuous finite element spaces where discontinuous Galerkin/finite volume methods can be 
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used. We shall use the flux 1-form representation of velocity u -idS in Ap(f2) (to evaluate di- 
vergence) and the circulation 1-form representation of velocity v-dxn (to evaluate vorticity) 
in A2(fi). 



Projection operators. Following [HipOlal . iHipOlbj . we define Hodge star projections II from 
A d (Q) k to A2- fc (fi) given by 



/ 7 a *Uuj 
Jn 



f 

Jn 



7 ACJ, V 7 6Aj(fi), 



and require that the dual spaces are chosen such that II is invertible. This requirement 
somewhat limits the choice of spaces, the main options arising from finite element spaces 



that have the same structure as the primal/dual C-grid staggerings, for example Chr08 



suggests the hexagonal P1-RT0-P0 spaces for the primal mesh, and the triangular P1-N0-P0 
spaces for the dual mesh, which would have the same structure as the hexagonal C-grid. 

Primal and dual 5 operator. The operator 5 provides weak approximations to V and V 1 - in the 
primal space, as described in the previous section for the primal finite element formulation, 
and weak approximations to V 1 and V- in the dual space. We shall now see that the key to 
the formulation is that we can define a simple relationship via the projection IT between d 
in the primal space and 5 in the dual space, and vice versa. 

Mapping from d to 5. For u e Aj(fi), we have 

/ 7A*ndc<j = / 7Adw 
Jn Jn 

= - / (I7AW 
Jn 

= — I d 7 a *Uui 
Jn 

= f 7A*<fflw, V 7 eA^(n), 
Jn 

where we may integrate by parts in the second step since d7 and du are both well-defined. 
Hence, Tlduo = 5Uu, and the following diagram commutes: 



A 2 d (fi) 



11 



Aj(n) 



AS(n) 



Ag(n) 



For example, this means that we may start with u -idS e A^(fi), and either obtain the primal 
vorticity ( p 6 A°(f2) by directly applying 5, or by inverting Pi to get v-dx^ e A^(f2), applying 
d to get the dual vorticity QdS 6 A^(f2), and then projecting back to A°(fi) with n, i.e., 



S(u JdS) = ndn- 1 (n jdS). 
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Prognostic equations. Within this framework, we retain the same equation for D on the 
primal grid, and modify the Coriolis term in the u equation as follows: 



fujd5+ f wjdSAUQ-dx n - f d(wJdS)A*(g(D-b)+K)dS = 0, V™ j dS e A^(fi), 
dt Jn Jn 



— f(/>dSA*(Dds)- [d*(<j)dS)/\u-\DdS+ [ *((/>dS)/\u-\DdS = 0, V^dSe Aj(ft), 

d t Je J e -/9e 

with Q • da?n ^ Ai(O). In our shorthand notation, we can rewrite the velocity equation as 

f) 

(12) 



^-u j dS + UQ ■ dx n + 5 (g(D - b) + K) = 0. 



Primal vorticity. As before, the primal grid vorticity ( p is defined by = 5u J d 5, and 
application of 5 to Equation (fl2"]l gives 



■^-cfri, jdS + OTQ-dxn = 0, 

at 



since £ 2 = 0, and hence 



j? 
94 



— C P + 5nQ-da; n = 0. 



(13) 



Dual vorticity. We now introduce an equivalent dual grid vorticity (ddS e Ag given by 
IlC^dS = ( p . Substituting u j dS = Tlv ■ dx^ with v ■ dxn 6 A^(f2), we obtain 

nC«jd5 , = Cp = <y«-id5 = <mt;-dajn = nd(t;-da;n) > 

and hence Q = d(v ■ dx^) by invertibility of II. Substitution into Equation ( 1T31) and applica- 
tion of the commutation relations for II gives 





and hence 



— nCd + ndQ-da; n = 0, 

at 



d 

—(d + dQ-dx n = 0, 
at 



by invertibility of II. It remains to determine a suitable Q. 

Dual potential vorticity. This enables us to define a potential vorticity q with qdS e A\(Q), 
given by 

f (f)dS A*(qDdS) = f <j)dS A*(( d + f)dS, \/(f)dSe k 2 d (Q). 
Jn Jn 

We then write down the potential vorticity equation on dual cell e^: 

^- f (pqDdS- f d(pAFjqdS+ f (/)AqFJdS = 0, V^eA^fi), 

dt -/e d Je d Jde d 

noting that if q is constant then we obtain 



r (f) ^ Dd s=- f 4> q 

Je d at Je d 



A 



\ 



|-Dd5 + dFjdS 

dt 



e A2(n), 
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and hence 

and we have mass consistency. 



Dual vorticity flux. Identical to the mass flux reconstruction on the primal grid, we can then 
find Q e Aj(fi) such that 

%-£ d dS + d(Q-dx n ) = 0, 
at 

and hence 

dC 3 

-^ = -U( d dS = -Ud(Q-dx n ) = -5U(Q-dx n ) 

which means that 

2-8u JdS + SU(Q-dx n ) = 0, 

at 

which is consistent with the primal velocity equation 

^-u j dS + Il(Q- dx n ) -5(gD + K) dS = 0. 
at 



5. Summary and outlook 

In this paper, we used the finite element exterior calculus framework to develop two for- 
mulations for the shallow-water equations, a primal formulation that is defined on a single 
mesh where divergence is defined in the strong form but vorticity must be evaluated weakly 
using integration by parts, and a primal-dual formulation that makes use of a dual mesh 
where vorticity can also be computed in the strong form. Both of these formulations ad- 
dress the list of desirable properties given in the introduction, and they additionally have a 
conserved diagnostic potential vorticity. Both formulations provide a way to control oscilla- 
tions in the divergence-free component of the velocity field (the component that dominates in 
large scale balanced flow in the atmosphere) by ensuring that the potential vorticity remains 
mass-consistent and oscillation-free. In the primal-dual case this can be achieved since the 
potential vorticity is diagnosed on a discontinuous space where discontinuous Galerkin/finite 
volume methods can be used to provide stable shape-preserving potential vorticity fluxes. 
In the primal case, the potential vorticity is computed in a continuous finite element space, 
but streamline-upwind Petrov-Galerkin methods with discontinuity capturing are compatible 
with the framework and can be used to provide stable potential vorticity fluxes with minimal 
oscillations. 

This work is part of the UK GungHo Dynamical Core project, which is a NERC/STFC 
collaboration between UK academics and the UK Met Office to design a dynamical core 
for the Unified Model that will perform well on the next generation of massively parallel 
supercomputers. In Phase 1 of the project, one of the main goals is to determine the horizontal 
discretisation that will be used, with the shallow-water equations on the sphere providing an 
environment to investigate this. The aim is to develop discretisations on a pseudouniform 
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grief] that have all of the desirable properties listed in the introduction, whilst maintaining the 
accuracy of the current model. Numerical accuracy is crucial since it reduces grid imprinting 
(structure in the numerical errors that reflects the structure of the grid, e.g. larger errors 
near the corners of a cubed sphere). This work o pens u p a number of possibilities that could 



be sufficiently accurate for operational use. In |CS12| it was shown that to avoid spurious 
mode branches it is necessary to select finite element spaces that have a 2:1 ratio of velocity 
DOFs to pressure DOFs, which suggests the BDFM1 space on triangles with an icosahedral 
mesh in the primal formulation or RTO on quadrilaterals with a cubed sphere mesh in the 
primal-dual formulation. There is an argument to be made that spurious Rossby mode 
branches arising from increasing velocity DOFs relative to this ratio are not harmful since 
they have very low frequencies and will just be passively advected by the flow. This suggests 
the BDM1 space on triangles in the primal formulation or the RTO space on hexagons in 
the primal-dual formulation, which both have a 3:1 ratio. The next steps in this work are 
to analyse the numerical convergence and dispersion relations of all of these schemes and to 
benchmark them against the usual suites of testcases and against solutions from the Unified 
Model formulation. 
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a primal-dual finite element formulation, Darryl Holm for useful guidance and comments on 
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